'''
Copyright: 
Descripttion: 
version: 
Author: chengx
Date: 2020-12-23 14:43:44
LastEditors: chengx
LastEditTime: 2021-05-29 15:19:08
'''
'''
本代码是一种高光谱图像波段选择算法
采用ICA(独立分量分析)方法。
Reference: "Band Selection Using Independent Component Analysis for Hyperspectral
Image Processing", Hongtao Du, et al. 2003
Inputs:
        Should be given:
            A hyperspectral data-set in ENVI format
        Will be asked via the console:
            Number of component
            Number of the bands required to be selected
The number of component is determined by the user
It can be evaluated to minimise the difference
between the original band set (X) and multiplication of mixing matrix (A_) and the source (S_)



可以通过“assert”函数和改变“atol”和“rtol”参数(第62行)来进行评估，以使原始带集(X)和混合矩阵(A_)与源的乘法(S_)之间的差异最小。
通过‘assert’测试函数的参数值越小，混合矩阵的估计就越好
'''

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from sklearn.decomposition import FastICA

def myIcaBs(data):
    """
    Parameters: data :输入高光谱图像，shape为(长*宽,波段数)
    Description: 波段选择
    Returns: bandInd 得分较高的波段索引
    """

    X = data

    n_c = 30
    n_b = 15

    # FastICA需要对数据进行白化处理
    ica=FastICA(n_components=int(n_c), whiten=True) 
    S_ = ica.fit_transform(X)        # Reconstruct signals
    A_=ica.mixing_                 # 获得估计的混合矩阵

    # 计算非方阵矩阵的秩
    if A_.shape[1] != np.linalg.matrix_rank(A_):
        print('A does not have left inverse ')
    else:
        # A_的伪逆运算
        W=np.linalg.pinv(A_)    # W. transpose(X)=transpose(S_)
        assert np.allclose(A_, np.dot(A_, np.dot(W, A_)))  # to check the pseudo-inverse matrix
    B_W=np.sum(np.absolute(W),axis=0)   # 计算每个波段的权重
    sortB_W =np.argsort(B_W)   # 提取索引
    
    bandInd = sortB_W[-int(n_b):]
    return bandInd
